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ABSTRACT 

Mrk 231 is a nearby ultra-luminous IR galaxy exhibiting a kpc-scale, multi-phase AGN-driven outflow. This galaxy represents the best 
target to investigate in detail the morphology and energetics of powerful outflows, as well as their still poorly-understood expansion 
mechanism and impact on the host galaxy. In this work, we present the best sensitivity and angular resolution maps of the molecular 
disk and outflow of Mrk 231, as traced by CO(2-l) and (3-2) observations obtained with the IRAM/PdBI. In addition, we analyze 
archival deep Chandra and NuSTAR X-ray observations. We use this unprecedented combination of multi-wavelength datasets to 
constrain the physical properties of both the molecular disk and outflow, the presence of a highly-ionized ultra-fast nuclear wind, 
and their connection. The molecular CO(2-l) outflow has a size of ~ 1 kpc, and extends in all directions around the nucleus, being 
more prominent along the south-west to north-east direction, suggesting a wide-angle biconical geometry. The maximum projected 
velocity of the outflow is nearly constant out to ~ 1 kpc, thus implying that the density of the outflowing material must decrease 
from the nucleus outwards as - r -2 . This suggests that either a large part of the gas leaves the flow during its expansion or that 
the bulk of the outflow has not yet reached out to ~ 1 kpc, thus implying a limit on its age of ~ 1 Myr. Mapping the mass and 
energy rates of the molecular outflow yields M 0F = [500 - 1000] M© yr _1 and E kin 0 F = [7 - 10] x 10 43 erg s -1 . The total kinetic 
energy of the outflow is E kin 0 F is of the same order of the total energy of the molecular disk, E disk . Remarkably, our analysis of the 
X-ray data reveals a nuclear ultra-fast outflow (UFO) with velocity -20000 km s -1 , M UF0 = [0.3 - 2.1] M 0 yr _1 , and momentum 
load PuFo/Prad = [0.2 - 1.6]. We find E kin?UF0 ~ E kin 0 F as predicted for outflows undergoing an energy conserving expansion. This 
suggests that most of the UFO kinetic energy is transferred to mechanical energy of the kpc-scale outflow, strongly supporting that 
the energy released during accretion of matter onto super-massive black holes is the ultimate driver of giant massive outflows. The 
momentum flux P 0 f derived for the large scale outflows in Mrk 231 enables us to estimate a momentum boost Pof/Pufo - [30 - 60]. 
The ratios E k m,uFo/L bo i,AGN = [1 - 5]% and E kin?0 F/L bo i,AGN = [1 - 3]% agree with the requirements of the most popular models of 
AGN feedback. 

Key words. Galaxies: individual: MRK231 - Galaxies: active - Galaxies: evolution - Galaxies: ISM - Galaxies: kinematics and 
dynamics - Galaxies: quasars: general 


1. Introduction 

Luminous AGN are likely to play a key role in the evolution of 
their host galaxies through powerful winds and outflows. In par¬ 
ticular, during the earliest, dust obscured phase of AGN activity, 
characterized by a large amount of molecular gas available to 
fuel black hole accretion, the huge nuclear energy output enables 
the production of AGN-driven outflows. These can perturb and 
possibly expel most of the gas out of the host galaxy, therefore 
leading to the quenching of star-formation. Current theoretical 
models of galaxy evolution include AGN feedback to account 
for galaxy colours and the lack of a large population of local 

* Based on observations carried out with the IRAM Plateau de Bure 
Interferometer. IRAM is supported by INSU/CNRS (France), MPG 
(Germany) and IGN (Spain), and with Chandra and NuSTAR obser¬ 
vatories. 


massive star forming galaxies (SFG, Granato et al. 2001, 2004, 
Di Matteo et al. 2005, Menci et al. 2006, 2008, Bower et al. 
2006, Croton et al. 2006). 

Observations across the electromagnetic spectrum reveal that 
AGN winds are widespread, suggesting that they cover large 
solid angles and an extremely broad range of gas ionization pa¬ 
rameter (Elvis 2000 and references therein). In X-rays mildly 
ionized warm absorbers with velocities of -500-1000 km s -1 
are observed in more than half of AGN (Piconcelli et al. 2005; 
Blustin et al. 2005; McKernan et al. 2007). Furthermore, there is 
growing evidence for the existence of ultra-fast outflows (UFOs) 
via the detection of highly ionized Fe K-shell absorption arising 
from a gas moving at velocities up to 0.2-0.4c (Pounds et al. 
2003; Reeves et al. 2009; Tombesi et al. 2011, 2015 and refer¬ 
ences therein). Optical/UV Broad Absorption Line (BAL) sys¬ 
tems associated to ionized outflowing (-2000 -f -20000 km s -1 ) 
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gas has been observed in 20-40% of QSOs (Dai et al. 2008), with 
some remarkable examples reaching up to hundred of pc away 
from the nucleus (e.g., Borguet et al. 2013). Broad [OIII] emis¬ 
sion lines with velocities up to a few thousands km s -1 are com¬ 
monly found in AGN (Zhang et al. 2011; Shen & Ho, 2014), and 
ionized and neutral gas with similar velocities is found on galac¬ 
tic scale (1-30 kpc) in a few tens AGN from the local Universe 
up to z=2-3 (Cano-Diaz et al. 2012, Rupke & Veilleux 2011, 
2013; Harrison et al. 2014; Nesvadba et al. 2008; Harrison et al. 
2012; Genzel et al. 2014, Cresci et al. 2015). Radio observations 
have also detected many neutral atomic gas outflows with veloc¬ 
ities up to 1500 km s -1 (Morganti et al. 2005; Teng et al. 2013). 
Finally, fast (up to 1000-2000 km s -1 ) molecular gas, extending 
on kpc scales is found in a few dozen, mostly highly obscured, 
AGN in dusty SFGs (Feruglio et al. 2010, Cicone et al. 2014, and 
Veilleux et al. 2013 for compilations of molecular outflows). 

Two main scenarios for quasar-mode winds and feedback 
have been developed so far. The first is a two-stage mechanism 
(Lapi et al. 2005; Mend et al. 2008; Zubovas & King 2012; 
Faucher-Giguere & Quataert 2012), where a (semi)relativistic 
wind (such as the UFOs detected in the X-ray band), accel¬ 
erated by the AGN radiation, shocks against the galaxy ISM. 
The shocked gas can reach temperatures of 10 7 K or larger. Hot 
gas extended on kpc scales has been detected in galaxies with 
galactic-scale outflows (Feruglio et al. 2013, Wang et al. 2014, 
Veilleux et al. 2014). If the cooling time of the post shock gas 
is smaller than the dynamical time, t coo i < tdyn, then the gas 
can cool to - 10 4 K, emitting typical warm ionized gas lines 
([OIII], H a, etc). Further cooling leads to the formation of H 2 , 
which can further enhance the cooling rate, promptly leading to 
the formation of more complex molecules, and a multi-phase, 
galaxy-scale outflow with a momentum boost Pof > 10 x P ra d 
should develop. Simulations of AGN feedback triggered by cold 
precipitation produce indeed a multi-phase structure of the gas 
in the galaxy cores (Gaspari et al. 2013). 

In the second scenario, outflows are radiatively accelerated 
by much slower, ’gentle’ winds, in which the momentum boost 
is at most Pof < 5 x P ra d (Thompson et al. 2014). This may pre¬ 
vent gas clouds from escaping the halo gravitational potential, 
hence making them rain back onto the galaxy. Self-shielding of 
the galactic disk dusty clouds may further minimize AGN feed¬ 
back. 

The geometry of quasar-driven outflows as well as their en¬ 
ergy and momentum transport mechanisms are still poorly con¬ 
strained by observations. This prevents both a clean discrimi¬ 
nation between the two scenarios, and the development of de¬ 
tailed multi phase feedback models. This work aims at filling 
this gap by means of high spatial resolution CO mapping and 
deep X-ray spectroscopy of Mrk 231. Mrk 231 has been known 
for long as the closest ULIRG/QSO, with a total infrared lumi¬ 
nosity of 3.6 x 10 12 L 0 (Sanders et al. 2003), AGN bolomet- 
ric luminosity of 5 x 10 45 erg s -1 , and a star formation rate of 
140 M© yr -1 (Veilleux et al. 2009). It is an iron low-ionization 
broad absorption-line quasars (FeLoBAL), and it displays nu¬ 
clear neutral and ionized outflows in several optical and UV 
tracers (Lipari et al. 2009, Veilleux et al. 2014), and a powerful 
neutral outflow extended on at least 3 kpc (Rupke & Veilleux 
2011). This is tracing a wide-angle conical outflow proceed¬ 
ing at an uniform speed of -1000 km s -1 along the minor axis 
of the molecular disk. The first prototypical massive molecular 
outflow, extended on kpc scale, was discovered in this galaxy 
(Feruglio et al. 2010, Fischer et al. 2010, Aalto et al. 2012, 2014, 
Gonzalez-Alfonso et al. 2014). These observations leave several 
questions open, i.e. (i) what are the size and the geometry of the 


molecular outflow, (ii) whether and how this outflow is associ¬ 
ated with the atomic neutral and ionized ones, (iii) which is the 
momentum boost of the molecular outflow, (iv) is there evidence 
of any fast highly ionized nuclear wind in Mrk 231? 

This work is a study of the multi-phase winds in Mrk 231 
(z=0.04217, D l = 187.6 Mpc, scale 0.837 kpc/"). Specifically, 
Section 2 presents the (sub)millimeter observations of molecular 
gas traced by CO(2-l) and (3-2), obtained with the IRAM/PdBI. 
Section 3 presents an analysis of the kinematics of the nuclear 
region of Mrk 231, including the molecular disk and the out¬ 
flow, and spatially resolved mass flow rates and outflow kinetic 
power. Section 4 presents an analysis of archival X-ray data from 
Chandra and NuSTAR , that led us to detect a nuclear, highly 
ionized UFO. In Section 5 the relation between the spatially- 
extended molecular outflow and the UFO is discussed in the 
framework of models for energy-conserving large-scale winds 
driven by AGN. 


2. (Sub)-Millimeter Observations 

We use three data sets obtained with the IRAM Plateau de Bure 
Interferometer (PdBI). The first dataset (DS1) consists of ob¬ 
servations carried out with the most extended (A) 6 antenna ar¬ 
ray configuration, which has a maximum baseline length of 760 
m, at a frequency of 221.21 GHz, targeting the CO(2-l) line. 
Absolute flux calibration is based on the quasars 3C84 (13.2 Jy) 
and 1150+497 (2.6 Jy), which also served as amplitude/phase 
calibrator. The accuracy of the absolute flux scale is of the or¬ 
der of 20%, while the relative flux calibration between the 4 
available observations is better than 10%. The on source time 
of DS1 after calibration is 5.2 hours. The synthesized beam is 
0.48 by 0.36 arcsec, PA 73 deg (-400 pc at the redshift of the 
source z= 0.04217). The achieved visibility (thermal) noise is 
1.3 mJy/beam in 20 MHz channels (corresponding to 27 km 
s -1 ). Natural weighting and a simple cleaning algorithm {hog- 
bom) were used in order to minimize secondary lobes. The visi¬ 
bility tables and data cubes have been binned to a 20 MHz spec¬ 
tral resolution, which is our preferred trade-off between spectral 
sampling and sensitivity. 

The second data set (DS2) has been obtained by merging the 
visibilities of DS1 with previous observations with short base¬ 
lines reported in Cicone et al. (2012). The relative flux calibra¬ 
tion of the short and long-baseline data is based on the quasar 
3C84, which has varied from 7.7 Jy in 2010 to 13.2 Jy in 2013. 
Based on the light curve of 3C84 around 220 GHz we can expect 
that the two data sets have consistent flux scales at the 20% level. 
DS2 has an on source time of 8.7 hours, and l<x sensitivity 0.8 
mJy/beam in 20 MHz channels. Maps for DS2 were produced 
as detailed above, resulting in a clean beam of 0.9 by 0.8 arcsec, 
PA 57. 

The third data set (DS3) consists of an observation of CO(3- 
2) at 331.8 GHz done in January 2012, using the compact (D) 
array configuration. The maximum phase rms during the obser¬ 
vation was about 50 deg, water vapor around 1 mm, and sys¬ 
tem temperatures between 300 and 500 K. The quasar 1150+497 
(0.96 Jy at 331.9 GHz) was used for phase/amplitude calibration. 
MWC349 was used for flux calibration and provides an accu¬ 
racy of about 20% at this frequency. The synthesized beam is 
1.5 by 1.3 arcsec, with a PA of 84 deg. The lcr sensitivity is 8.0 
mJy/beam per 20 MHz. 

DS 1 has the best angular resolution and will be used to locate 
and characterize the kinematics of the inner - 1 arcsec of the 
system. DS2 has the best sensitivity and will be used to map 
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Table 1 . Plateau de Bure Interferometer observations of Mrk 231 used in this work. 


Name 

Date 

Array 

Freq. 

[GHz] 

Line 

On source time 
[hr] 

Clean beam 
[arcsec] 

Sensitivity 

[mJy/beam/20MHz] 

DS1 

Feb 2013 

A 

221.21 

CO(2-l) 

5.2 

0.48 x 0.36 

1.3 

DS2 

Sep-Nov 2010/Feb 2013 

ACD 

221.21 

CO(2-l) 

8.7 

0.9 x 0.8 

0.8 

DS3 

Jan 2012 

D 

331.8 

CO(3-2) 

2.0 

1.5 x 1.3 

8.0 



-1000 -5Q0 0 500 1000 1500 2000 

L5R velocity (km/s) 

Fig. 1 . Continuum-subtracted spectra of CO(2-l) from DS1 (top panel) and of CO(3-2) (bottom panel). The plots display the full 
velocity range achieved by the observed bandwidth. Zoomed views are shown in the insets. 


fainter and more diffuse emission. Table 1 summarizes the main 
parameters of the observations. 

2.1. 1.4 mm data 

The 1.4 mm continuum was estimated by averaging the visibil¬ 
ities in two spectral regions free of lines in the range -3000 to 
-2000 km s"\ and 1500 to 1600 km s _1 . We tested different 
channel ranges to average the continuum visibilities, mainly on 
the blue side (1500 to 3000 km/s negative velocities). The fea¬ 
tures seen at velocity -600 out to -1000 km/s, and those out 
to +1000 km/s, persist using different channel ranges to esti¬ 
mate the continuum. The significance of spectral features be¬ 
yond these limits are more difficult to assess. In particular, on 
the red side, where the spectral range is limited to +1800 km/s, 
the estimate of the continuum may be more uncertain, and fea¬ 
tures redder than +1000 km/s could be spurious. A point source 
fit to the visibilities averaged as explained above, gives a 1.4 mm 
continuum flux density of 43.8 + 0.5 mJy (Table [2]), in agree¬ 
ment both with the recent interferometric measurement of Aalto 
et al. (2014), done with similar baselines, and with Downes & 
Solomon (1998). 


The continuum visibility table, derived as explained above, 
has been subtracted from the total uv data to produce the line 
table and the corresponding map. Figure 1, upper panel, shows 
the CO(2-l) spectrum, extracted in a circular region of diameter 
1.5 arcsec centered on the emission peak on the map. The peak 
flux of CO(2-l) is 0.73 Jy, from a Gaussian fit of the spectral 
line, and agrees within the flux accuracy (+%20) with previous 
interferometric measurements of Cicone et al. (2012), Downes 
& Solomon (1998) (which had slightly larger beam), and with 
the flux measured with the IRAM 30 m dish (Solomon et al. 
1997). This suggests that little emission is filtered out by long 
baselines and that most of the CO emission occurs within a 1.5 
arcsec central region. The line profile is similar to that found in 
previous CO(l-O) and (2-1) observations (Feruglio et al. 2010, 
Cicone et al. 2012). CO(2-l) emission is detected out to about 
+ 1000 km s -1 . 

Figure 2 (upper panels) shows the channel contour maps of 
CO(2-l) emission integrated around the systemic velocity (+300 
km s _1 ), in the blue (-400 to -1000 km s _1 ), and red (400 to 1000 
km s -1 ) sides of the line. Both the systemic and the red/blue 
shifted emission are more extended than the synthesized beam 
and have a size up to ~ 1.5 arcsec. We performed visibility fitting 
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Fig. 2. Channel maps of CO(2-l) (DS1, upper panels) and CO(3- 
2) (lower panels) in three different velocity ranges (indicated 
by the labels). The contours start from 4 cr and are in steps of 
lcr for the blue and red wing maps (right and left panels for 
both CO transitions, respectively). Negative contours are shown 
by dotted lines, and are in steps of lcr, starting at -lcr. In 
the systemic component map (middle panels), contours are in 
steps of 10cr starting from 10<x, for the sake of clarity. The syn¬ 
thesized beams are shown by hatched areas for each data set. 
lcr = 0.14 Jy km s' 1 for CO(2-l), = 0.72 Jy km s" 1 for CO(3-2). 
The cross marks the position of the 1.4 GHz continuum as mea¬ 
sured by the VLBI, a =12:56:14.2339, 8 =56:52:25.237 (J2000, 
Ma et al. 1998). 


shifted emissions are instead both offset from the AGN position 
in the south-western direction by (-0.2, -0.3) and (-0.2, -0.1) 
arcsec, respectively. These offsets are about ~10 larger than the 
position errors of our data, i.e., 0.02" for the blue and 0.01" 
for the red component, respectively (Table [ 3 ]). The positions of 
the wing emission are consistent with those seen in HCN(1- 
0) (Aalto et al. 2012), and with those previously measured by 
Cicone et al. (2012^] 

We discuss in the following the possibility that the spectral 
feature seen between -1000 and -700 km s -1 traces emission 
from a species different from CO. Specifically, this spectral fea¬ 
ture peaks at a frequency of 221.864 (-860 km s -1 ), which cor¬ 
responds to that expected for 13 CS(5 - 4) (rest frame frequency 
231.220 GHz). 13 CS(5-4) was first detected in an extragalactic 
source in the 1 mm survey of Arp220 (Martin et al. 2011), and 
traces very dense gas. Its critical density is, in fact, of the order of 
10 6 cm -3 . The integrated emission is detected with an 8<x signif¬ 
icance. Its FWHM is about 150 km s -1 , and the zero spacing flux 
is 2.5 ±0.3 mJy, by fitting an unresolved source. 12 CS(5-4) is un¬ 
detected in the IRAM 30 m spectrum down to a lcr sensitivity of 
~1 mK (8.7 mJy, Zhang, PhD thesis). The isotopic ratio 12 C/ 13 C 
in Mrk 231 is about 40, based interferometric data of 12 CO and 
13 CO J=l-0 (Feruglio et al. in preparation), which agrees with 
the lower limit given by Glenn & Hunter (2001). Assuming that 
13 CO comes from approximately the same region as 12 CO(1-0), 
we can apply this isotopic ratio to derive the expected strength of 
13 CS(5-4), which would be < 0.2 mJy. This suggests that the fea¬ 
ture can hardly be identified with 13 CS(5-4), and we conclude, 
therefore, that it is due to approaching CO(2-l). 


of these three CO(2-l) line components, systemic, blue and red- 
shifted emission, after having averaged each visibility set in the 
corresponding spectral range (similarly to Cicone et al. (2012), 
and references therein). The positions, fluxes and sizes so de¬ 
rived are reported in Table [3] for each component, together with 
the respective spectral range and the adopted source model. The 
visibility fits are independent of the synthesized beam and of the 
CLEAN algorithm, and there is no need to deconvolve the appar¬ 
ent size on the map. The integrated fluxes agree within the errors 
with the measurements of Cicone et al. (2012). The FWHM size 
of the systemic gas component is 0.86 + 0.01 arcsec, correspond¬ 
ing to a physical scale of ~ 700 pc. The bulk of the blue and red 
wing emission of CO(2-l) is compact, having FWHM sizes of 
about a beam by fitting a with Gaussian function. A fainter spa¬ 
tially extended structure is seen in the channel maps of the red 
and blue wings (Fig. 3), which can hardly be fitted by a simple 
gaussian model, and, in the case of the receding component, is 
extended out to at least 1"(~ 1 kpc) north-east off the AGN. This 
component is also seen in HCN(l-O), on similar spatial scales 
and orientation (Aalto et al. 2012). The positions of both the 
systemic CO emission and of the 1.4 mm continuum are consis¬ 
tent with that of the 1.4 GHz radio continuum measured by the 
VLBI, which traces the AGN emission. The blue and the red- 

Table 2. 1.4 and 0.9 mm continua visibility fit parameters 


Wavelength 

Model 

RA 

DEC 

S v,zero—s pacing 

1.4 mm 

point source 

12:56:14.23 

56:52:25.24 

43.8 ± 0.5 mJy 

0.9 mm 

point source 

12:56:14.24 

56:52:25.23 

72.0 + 2.0 mJy 


Note. - Errors: flux errors are statistical, additional +20% should be accounted for flux cali¬ 
bration. 


2.2. 0.9 mm data 

The 0.9 mm continuum has been estimated by fitting in the uv 
plane visibilities averaged in spectral windows. We have tested 
the continuum estimate and subtraction by using three differ¬ 
ent spectral windows for visibility averaging. In particular, av¬ 
eraging the spectral windows -1180 to -1000 km/s, plus 1400 
to 2000 km/s produces over-subtraction of the continuum red- 
wards of the line. We find that the best overall continuum sub¬ 
traction results from using the spectral window 1400 to 2000 
km/s, which gives a continuum flux density of 72.0 ± 2.0 mJy 
(Table [2]). This is our adopted best solution. The continuum vis¬ 
ibility table, produced as explained above, has been subtracted 
from the total visibilities in order to produce the line emission 
visibility table, which has been used to map the CO(3-2) emis¬ 
sion (Figure 2, bottom panels). 

The CO(3-2) spectrum, extracted from a circular aperture of 
3 arcsec diameter around the map peak, is shown in Figure 1 
(bottom panel). We measure a flux density of 658 ± 5 Jy km s -1 
for the CO(3-2) systemic component. The single dish flux from 
Wilson et al. (2008), 480+100 Jy km s -1 , is slightly smaller than 
our estimate. A high velocity component is detected in CO(3- 
2) out to speeds of +1000 km s _1 , similarly to CO(l-O) and 
(2-1). The CO(3-2) profile is similar to that of CO(2-l), (1-0), 
showing an enhanced red wing compared to the blue one. We 
acknowledge that the systematic uncertainty on the continuum 
measurement can affect the significance of the faint features seen 
at large velocities (in particular blue-wards of the line at veloc¬ 
ity < -1000 km/s), but features within +800 km/s from the line 
peak are persistent by adopting any of our tested spectral win- 


1 Note that in Cicone et al. (2012) the phase reference center was 
wrongly quoted as the VLBI position, which produces an apparent dis¬ 
crepancy in the location of the outflow between that and this work. 
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Fig. 3. Channel maps of CO(2-l) from DS1. Contours are in steps of 5<x, starting at 5cr. lcr = 0.09 Jy km s 1 . 


dows. The positions of the systemic, the red and blue-shifted 
components of CO(3-2) are consistent with those derived for 
CO(2-l), within the uncertainty associated to the larger beam 
and to the coarser uv coverage of DS3 (Table [3]). 

In the following we discuss the possible contribution of 
H 13 CN(4 - 3) emission to the red wing of CO(3-2). Sakamoto et 
al. (2009) found a ratio CO(3-2)/H 13 CN(4-3)= 10inNGC4418, 
which implies severe contamination of the CO(3-2) red wing for 
that source. A ratio of ~ 10 is excluded for Mrk 231 by the ob¬ 
served spectrum already (Fig. 1). Costagliola et al. (2011) mea¬ 
sured for Mrk 231 a H 13 CN(1 - 0) = 1.7 Jy km s -1 at a 3<x 
level. From their spectrum, the line appears somewhat broader 
than the HCN(l-O), suggesting that it can be blended with some 
other nearby lines. Indeed, emission from SO (86.1 GHz), and 
S0 2 (on the red side at 86.639 GHz) are expected within ±300 
MHz . These transitions are as bright as H 13 CN in the Orion and 
SgrB2 star forming regions (Turner et al. 1989), which can in 
principle also be the case in Mrk 231. Therefore the measure of 
Costagliola et al. (2011) must be considered as an upper limit. In 
fact, this would provide a ratio HCN/H 13 CN(1 -0) ~ 8, based on 
the strength of HCN(l-O) measured by Aalto et al. (2012), i.e., a 
factor of five smaller than the isotopic 12 C/ 13 C ratio of 40 found 
by us. Adopting an isotopic ratio of 8 and HCN(4 - 3) = 65 + 13 
Jy km s -1 (Papadopoulos et al. 2007), we get H 13 CN(4 - 3) = 8 
Jy km s -1 , meaning that this would contribute 25% of the inte¬ 
grated flux in the red wing of CO(3-2). Adopting instead a ra¬ 
tio 13 CO/ 12 CO = 40, the contribution to the CO(3-2) red wing 
would be less than 5%. 

The estimates given above assume that H 13 CN(4 - 3) is 
optically thick. Should this not be the case, its flux may be 
even smaller. Modeling with RADEX non-LTE radiative trans¬ 
port models (van der Tak et al. 2007) to reproduce flux of CO 
and HCN, suggests instead that H 13 CN(4 - 3) is optically thin. 
In fact, assuming Av = 60 km s -1 and typical temperatures 
and column densities (Tki n = 60 K, Nr = 1 x 10 16 cm -2 for 


HCN and 1 x 10 15 cm' 2 for H 13 CN), we get r ~ lx 10“ 2 for 
H 13 CN(4-3). In this case the flux of H 13 CN(4-3) would be very 
small and the contamination of the red wing of CO(3-2) negli¬ 
gible. The maps of the CO(3-2) give an observed temperature of 
T = 0.2 K in the red wing. The HCN(4-3) flux (65 Jy km s -1 , 
Papadoupulos et al. 2007), implies 0.325 Jy for a line width of 
200 km s -1 . Ignoring the frequency difference between 340 and 
330 GHz, than we would get 5.7 K/Jy for the PdBI data, thus 
T(HCN(4 - 3) = 1.85 K (i.e. 0.325 Jy). Modeling with RADEX, 
and assuming N(H 13 CN) ~ 1 x 10 15 cm -2 , n(H 2 ) ~ 1 x 10 4 cm -3 , 
and T^ = 60 K, we get T(H 13 CN(4 - 3)) ~ 2.5 x 10 2 K, and 
S(H 13 CN(4-3)) = 4.4 mJy, which in turn would give 0.88 Jy km 
s -1 integrated intensity for a line width of 200 km s -1 . In con¬ 
clusion, a strong contamination from H 13 CN(4 - 3) is unlikely 
and most of the detected flux in the red wing of CO(3-2) is from 
high speed CO(3-2). 

3. Results of (sub)-mm observations 

3.1. Molecular disk 

We built moment 1 and 2 maps of the CO(2-l) emission us¬ 
ing DS1. The Moment 1 (velocity) map of CO(2-l) in Figure |4j 
left panel, shows that rotation occurs in approximately the east- 
west direction, as previously measured by Downes & Solomon 
(1998), and the molecular disk extends out to a radius about 800 

Table 4. Fit parameters derived for the molecular disk. 


Component 

ro 

PA 

Line-of-sight Inclination 


[arcsec] 

[deg] 

[deg] 

Outer ring 

0.24 

-12 + 5 

36 + 10 

Inner disk 

0.12 

84 + 5 

58 + 10 
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Table 3. Visibility fits of the CO(2-l) and CO(3-2) 


Component 

Vel. range 
[km s -1 ] 

Dataset 

Model 

RA 

DEC 

S v,zero-spacing 

[mJy] 

JSydv 

[Jy km s 1 ] 

Source Size FWHM 
[arcsec] 

CO(2-l) Systemic 

-300++300 

DS1 

expon. disk 

12:56:14.230 [+0.003”] 

56:52:25.200 [+0.002”] 

311.0 + 0.6 

186.6+3.8 

0.86+0.01 

CO(2-l) Blue 

-400 +-1000 

DS1 

gauss 

12:56:14.22 [+0.02”] 

56:52:25.06 [+0.02”] 

6.9+0.1 

4.2+0.5 

0.35+0.07 

CO(2-l) Red 

400 + 1000 

DS1 

gauss 

12:56:14.23 [+0.01”] 

56:52:25.04 [+0.01”] 

12.3+0.1 

7.4+0.6 

0.44+0.05 

CO(3-2) Systemic 

-300++300 

DS3 

expon. disk 

12:56:14.24[±0.003”] 

56:52:25.19[±0.003”] 

1096 ± 9 

658 + 5.0 

1.150 + 0.015 

CO(3-2) Blue 

-400 +-1000 

DS3 

gauss 

12:56:14.24[±0.09”] 

56:52:25.00[±0.08”] 

26.0 ± 5.8 

15.6 + 3.5 

0.6 ± 0.4 

CO(3-2) Red 

400 + 1000 

DS3 

gauss 

12:56:14.23[±0.05”] 

56:52:25.22[±0.05”] 

51.6 + 6.0 

31.0 + 3.6 

0.8 ± 0.2 


Note. - Errors: flux errors are statistical, +20% should be accounted for flux calibration. The physical scale is 0.837 kpc/". 



RA RA 

Fig. 4. Moment 1 map (left panel, levels are 20 km s 1 each) and Moment 2 map (right panel, levels are 25 km s 1 each) of the 
systemic component of CO(2-l). Cross: AGN position from Ma et al. (1998). 


pc. The S-shaped pattern (Fig. [4]) indicates that the kinematics in 
the inner disk deviates from a simple east-west rotation, and sug¬ 
gests the presence of an inner warped disk on scales of < 0.2", 
tilted with respect to the main rotation plane. The Moment 2 map 
(Figure |4j right panel) is peaked close to the AGN, and shows 
branches of emission in the north, south and west directions, also 
indicating an additional kinematical component. 

By fitting an exponential disk model to the visibilities in the 
velocity range -300 to 300 km s -1 (systemic component) we de¬ 
rive a similar size of the molecular disk of 720 pc FWHM (Table 
4). The disk size derived from CO (3-2) is about 900 pc. To better 
constrain the properties of the molecular disk, we first fit the flux 
and velocity maps of the molecular disk (systemic component) 
using the kinematical disk model by Gnerucci et al. (2011a). We 
assumed that (i) the gas is circularly rotating in a thin disk, (ii) 
the gravitational potential depends only on the disk mass, and 
(iii) the disk surface mass density distribution is exponential, i.e. 

£(r) = S 0 e~ rlr ° (1) 

where r is the distance from the center and ro is the scale radius. 
Following Cresci et al. (2009), the adopted surface brightness 
profile for both components is exponential 

/(r) = I 0 e~ rlra (2) 

The disk model and fit results are shown in Fig. [5] top panel. 
The model is then convolved with the beam of the data. The 
flux and velocity maps are fitted with a minimizing-^ 2 method. 
Residuals due to the inner S-shaped feature are seen in the ve¬ 
locity map (Fig. [5] top right panel). 


In order to account for the inner tilted disk we fit the veloc¬ 
ity map with two components, i.e., an inner disk and an outer 
ring with the same exponential profile but different scale radii 
ro, and allowing for different inclinations (/) and position angles 
(PA). Position angles are measured in deg and are positive from 
east via north. This fit gives for the outer ring a scale radius of 
ro = 0.24 arcsec, a PA of the line of nodes of -12 deg, and incli¬ 
nation of i = 36+10 deg (including a dominant contribution from 
systematic uncertainties). For the inner disk, we find a scale ra¬ 
dius of r 0 = 0.12, a PA= 84 deg, inclination =58 + 10 deg (Table 
[4]). The size (radius 200 pc), inclination and orientation of the ve¬ 
locity gradient (nearly north-south) of the inner warped disk of 
CO(2-l) match well those seen in the vibrationally excited tran¬ 
sition HCN(3-2) V 2 = 1 (Aalto et al. 2014), and that seen in the 
OH megamaser observation of Klockner et al. (2003). 

The dynamical mass enclosed within 400 pc from the AGN 
is 1.6 ^q 5 x 10 9 M q (with the errors dominated by the uncertainty 
on the inclination). Downes & Solomon (1998) found 1.27 x 
10 10 M© for an inclination of 10 deg (based on the ratio of the 
two semi-axis of the disk, and, therefore, probably affected by 
even larger uncertainties), with our best fit inclination of 36 deg 
it would be a factor of 10 smaller and, thus, in agreement with 
our measurements. 


3.2. Molecular outflow 

Table [3] reports the results of the visibility fits of CO (2-1) and 
(3-2). The FWHM size of the high velocity wings of CO(2-l) is 
0.35" (0.3 kpc) and 0.44" (0.4 kpc) for the red and blue compo¬ 
nents, respectively. In addition we detect lower surface bright- 
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Fig. 5. From left to right: model, model folded with the beam, observed distribution and residuals. Top panels: model of the molecu¬ 
lar disk with an exponential disc function. Middle panels: model with two components, an outer ring and an inner exponential disk. 
Bottom panels: velocity dispersion of the two-component model. 



Fig. 6. Position-velocity plots with east-west (left panel, left to right), south-north (middle panel, left to right) and PA -45 deg, 
south-west to north-east (right panel, left to right) cuts, through the CO(2-l) peak from DS1. Contours levels are 3 to 15<x, l<x =1.3 
mJy/20 MHz. 
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Fig. 7. Position-velocity plots of CO(2-l) from DS2 along the same slices than in Fig. [6] Contours levels are 3 to 15<x, lcr =0.8 
mJy/20 MHz. 
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Fig. 8. Position-velocity plots of CO(3-2) along the same slices than in Fig. [b] Contours levels are 3 to 15<x, lcr =0.8 mJy/20 MHz. 



Fig. 9. Velocity maps of the blue (left) and redshifted CO emission (right). The velocity (in km s 1 with respect to the systemic 
velocity, i.e. the CO(2-l) peak, is shown in the colour bars. The cross marks the AGN VLBI position. 


ness extended emission due to receding (redshifted) gas out to at 
least 1 kpc north-east off the AGN (Fig. 3). The CO (3-2) out¬ 
flow component has a FWHM size of 0.6 (500 pc) and 0.8" (700 
pc) for the blue and red components, respectively. 

Position-velocity plots along three different directions are 
shown in Figures [6] to [8] For each data set, we extracted a hori¬ 
zontal (i.e. along the major axis of disk rotation, west-east cut), 
a vertical (south-north) and a -45 deg (south-west to north-east) 
cut. In all three data sets the rotation pattern of the outer disk 
is seen in the west-east cuts. Rotation is also visible along the 
other two directions in the high resolution data. Emission with 


velocity close to the systemic (±300 kms/s) is extended out to 
about 1 arcsec from the nucleus. 

The high velocity gas (v > 400 km s -1 ) is seen in all data 
sets, and along all the examined directions, suggesting that the 
velocity distribution is rather isotropic. However, both the high 
speed receding and approaching gas components extend south¬ 
ward of the nucleus, with the approaching one further elongated 
along the south-west direction (see middle and right panels of 
Figure 6). 

Figure [9] shows the Moment 1 maps of the red and blue 
CO(2-l) wings in the velocity ranges -500 to -1000 and 500 to 
900 km s -1 , respectively. In the blue component a velocity gradi- 



































C. Feruglio et al.: The multi-phase winds of Mrk 231 


ent is seen. The speed of the gas is-900 km s -1 at the position 

of the nucleus, and decreases to -600 km s -1 at 0.4" south-west 
of the nucleus. This is consistent with the position-velocity plot 
along the diagonal cut (Fig. [6j right panel). The bulk of the gas 
receding at velocities of 700-800 km s -1 is also located at ~ 0.3 
arcsec south-west of the AGN (Figure[9] right panel), and shows 
a milder velocity gradient along approximately the same axis as 
the approaching gas. Emission with velocity of about 700 to 800 
km s -1 is also seen at (+0.8,+0.8") offset from the phase center, 
along the same diagonal direction where the blueshifted velocity 
gradient is found. We cannot exclude that this feature is a resid¬ 
ual side lobe, but the absence of the corresponding side lobe in 
the opposite direction, suggests that this feature might be real. 

3.3. Mapping M 0 f and E kin ^ 0 F 

We fitted the pixel-by-pixel spectra with three Gaussian func¬ 
tions, in order to simultaneously account for both low veloc¬ 
ity (rotating disk) and high velocity gas (outflows). We define 
the maximum velocity, v max , as the shift between the veloc¬ 
ity peak of broad component and the systemic velocity plus 
2<r, where cr is the velocity dispersion of the broad gaussian 
component iy max = velocity shift broa d + 2cr broa d , Abroad = 
FWHM broa d/235). The v max map is show in Figure [T0[ Based 
on these three component fits, we attempted to map the mass 
outflow rate and kinetic energy. In previous works (e.g. Feruglio 
et al. 2010) the mass outflow rate has been computed using the 
continuity fluid equation 

Mof — Cl R of pof v max (3) 

where pof is the average mass density of the outflow, v max is 
the maximum velocity of the outflowing gas, and R of is the ra¬ 
dius at which the outflow rate is computed, and Cl is the solid 
angle subtended by the outflow. Assuming a spherical sector, 
Pof ~ 3M 0 f/QR 3 o F , then, M 0 f = 3 X Vmax x M 0 f/Rof- In 
this formula, Rof is the distance from the nucleus (RA,DEC= 
12:56:14.23, 56:52:25.20). Accordingly, M 0 f represents the in¬ 
stantaneous outflow rate of the material at the edge Rof (he., it 
is a local estimate) and it is three times larger than the total out¬ 
flow mass divided by the time required to push this mass through 
a spherical surface of radius Rof- This estimator does not de¬ 
pend on the solid angle Cl subtended by the outflow. Figure |TT] 
is the map of the outflow rate Mof derived by averaging the CO 
flux per beam of the broad gaussian components in quadrants of 
increasing radii (south-west, south-east, north-east, north-west 
off the AGN), and by adopting aco - 0-5 to convert from line 
luminosity to gas mass (Weiss et al. 2001). Mof in Figure [TT[ 
therefore, should not be considered a pixel by pixel estimator of 
the mass outflow rate, but rather the value in square sectors at 
increasing radii. 

The outflow originates from the nucleus, and Mof decreases 
in all directions with increasing distance from the nucleus. The 
outflow expands preferentially toward south-west, out to a dis¬ 
tance of 0.4 arcsec from the nucleus. 

3.4. M of , E k m,oF and M 0F I M disk radial profiles 

In the following we study the mass outflow rate, kinetic power 
and M 0 F/M d i sk as a function of the distance from the nucleus. 
Since the outflow is located near the nucleus and rather symmet¬ 
ric, for the sake of simplicity we adopted an one dimensional av¬ 
erage. We used DS1 to map the region at < 0.5 kpc with the best 
angular resolution, and DS2 to constrain any fainter high speed 



2.4e+02 3.1e+02 3.8e+02 4.5eF02 5.2e+02 5.9e-l-02 6.6e4-02 7.3e+02 8e+02 __ 

Fig. 10. Maps of v max = velocity shiftbroad + 2cr hroad (colour scale 
in km s -1 ). The pixel size is 0.1 arcsec. The cross marks the AGN 
VFBI position. 



475 650 827 1002 1178 1353 1528 1705 1880 _ 

Fig. 11. Map of Mof (the scale in M© yr -1 is reported at the 
panel bottom). Mof is computed by averaging the CO flux of 
the broad gaussian components in quadrants of increasing radii 
(south-west, south-east, north-east, north-west squared sectors 
off the AGN), and then converting to M(H2) using a conver¬ 
sion factor of 0.5. The clean beam is shown as an open ellipse. 
The physical scale is 0.1 arcsec/pixel. The cross marks the AGN 
VFBI position. 


emitting gas further out from the nucleus. We extracted spectra 
from box regions centered on the nucleus and with increasing 
sides. We then fitted the spectra with three Gaussian functions 
(or four when the systemic line needs two components) to model 
both the rotating disk and the outflow components, and com¬ 
puted the mass outflow rate as detailed above. The gas masses 
of the disks (traced by the systemic CO line) and outflow (broad 
CO components) have been calculated by integrating the best fit 
models and by using aco = 0.5. 

Figure [T2| compares the best fit CO(2-l) emission line pro¬ 
files integrated in a square region around the nucleus (using 
DS1) with the integrated profiles from two adjacent square an¬ 
nuli at increasing distances from the nucleus (using DS3). The 
rms is ~ 1.3 x 10 -3 Jy for all the spectra. For the receding gas, 
we detect high speed (800 km s -1 ) gas out to > 1 kpc, and a 
deficit of gas with intermediate velocity (300-500 km s -1 ) at 
>0.5 kpc. For the approaching gas, this deficit is only seen at 
>0.9 kpc. The outflow mass rate Mof, the v max , the kinetic en¬ 
ergy rate E^of = 0.5xM O FXv^ ax , and the ratio of outflow mass 
and molecular disk mass M 0 F/Mdisk are shown in Figure |~j~3] as 
a function of the distance from the nucleus (error-bars represent 
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Fig. 12. CO(2-l) emission line profiles extracted from square 
regions at different distances from the nucleus, as indicated by 
the colour-coded labels, and their multi-Gaussian best fit. 

the statistical errors only). Specifically, the histograms represent 
integral quantities out to a given radius, while the points rep¬ 
resent the local mass outflow rate in two annuli, computed by 
measuring the mass density and the outflow mass within the an¬ 
nuli. The integral mass outflow rate M 0 f is ~ 1000 M© yr -1 
within 400 pc from the nucleus, and 500-700 M© yr -1 out to 
~ 1 kpc. It is worth noting that the local mass outflow rate is 
about 500 M© yr -1 within ~ 800 pc, while it drops to a few 
tens M© yr -1 at > 1 kpc. The v max and the integral Eki n ,OF 
of the outflow remains nearly constant out to ~ 1 kpc, with 
Ekin,OF = 7 - 10 x 10 43 erg ~ 1-2% of the AGN bolometric lu¬ 
minosity (Figur4T3| middle panel, see Section 5 for a detailed 
discussion). Finally, Figure |T3) right panel, shows that the out¬ 
flow carries ~0.2-0.25 of the total disk mass out to ~ 1 kpc, 
while the outflow mass drops to less than 10% of the disk mass 
at > 1 kpc. 

4. X-ray observations 

4.1. X-ray data reduction 

During the last three years Mrk 231 has been target of new, 
sensitive X-ray observations. Specifically, Chandra observed the 
galaxy for 400 ks in August 2012 (Veilleux et al. 2014), while 
NuSTAR has observed it twice, in August 2012 and May 2013 
for a total of about 70 ks (Teng et al. 2014,T14 hereafter). These 
data have dramatically changed our understanding of the X-ray 
emission from Mrk 231. Previous broadband, non-focusing X- 
ray observations performed with BeppoSAX and Suzaku detected 
a ~3<x excess in the band above 10 keV which has been inter¬ 
preted as evidence of nuclear continuum emerging after trans¬ 
mission through a Compton thick absorber (Braito et al. 2004), 
most likely with a variable covering factor (Piconcelli et al. 
2013). This scenario has not been confirmed by the unprece¬ 
dented angular resolution ultra-hard (>10keV) X-ray NuSTAR 
observations presented by T14. They did not report any hard X- 
ray excess and revealed that Mrk 231 is therefore intrinsically 


X-ray weak, with a 2-10 keV luminosity of 4 x 10 42 erg s -1 . 
The best-fit model of the contemporaneous Chandra + NuSTAR 
spectrum consists of flat (r « 1.4) power-law continuum emis¬ 
sion modified by a patchy, Compton-thin absorber, plus a soft 
X-ray, starburst related, thermal emission. Furthermore, the deep 
Chandra observation has revealed the existence of a huge (~ 65 
x 50 kpc) soft X-ray halo around the central AGN which can 
be accounted for by two thermal emission components with kT 
~ 0.25 and 0.8 keV, respectively (Veilleux et al. 2014). Thanks 
to their high quality and sensitivity, these datasets also allow a 
detailed search for highly ionized fast or ultra-fast outflows seen 
in absorption against the nuclear X-ray emission. 

Chandra data were taken from the CXC archive. 
Specifically, we combined the 2012 long (400 ks, Observation 
ID 13947, 13948, 13949) observation with those performed in 
2000 and 2003 (153 ks in total, Observation ID 1031, 4028, 
4029, 4030, Gallagher et al. 2005). The combined data set has 
a total exposure time on source of 553 ks. Data were reduced 
using CIAO 4.5. We extracted a spectrum from a circular region 
of 3 pixel radius (1.5 arcsec) centered on the nucleus using 
the tool dmextract. A background spectrum was extracted 
from an annulus with inner and outer radii 1 and 2 arcmin, 
respectively. In extracting the background, the regions of the 
front-illuminated detector have been masked. We verified that 
varying the background extraction regions, however, has little 
impact, because the background counts are a small fraction of 
the source counts in the spectral region of interest (a factor of 
~ 1/500). Response matrices were computed using the tools 
mkwarf and mkrmf. The spectral analysis was performed in the 
0.5-10 keV energy range. Given the large number of counts 
and the requirement to use the x 1 statistics in our modeling, we 
binned the spectrum with a minimum of 40 counts/channel. 

The NuSTAR data were reduced with the pipeline 
NuSTARDAS version 0.11.1 and CALDB version 20130509 
with the standard settings (see T14 for details). The background 
counts are a factor of ~ 1/10 those of the source. Spectra were 
extracted for each observation and for the two NuSTAR tele¬ 
scopes FPMA and FPMB, using a circular region of 1 arcmin ra¬ 
dius. Spectra were binned with a minimum of 40 counts/channel 
as for the Chandra data set. Spectral bins between 3 and 79 keV 
were used in the fits, xspec 12.8.0 was used for the analysis. 

4.2. Discovery of a nuclear Ultra-Fast Wind 

We exploited these data sets to constrain any nuclear wind. 
Based on T14 results, we first fitted the Chandra spectrum and 
the four NuSTAR spectra with a model including Galactic ab¬ 
sorption along the line of sight, two thermal components, a 
power law component and a narrow emission line component, 
both reduced at low energies by photoelectric absorption (we 
used the xspec model zxipcf, i.e. a model including a partial cov¬ 
ering and ionized absorber). The total x 1 of this fit is 498.9 for 
468 deg of freedom (dof). Figure [14] shows the ratio between 
data and model. For plotting purposes the four NuSTAR spectra 
have been added together. Positive residuals are evident between 
6 and 6.5 keV (6.2-6.7 keV rest frame, consistent with neutral 
Fe K a emission) while negative residuals are seen around 7 keV. 
The feature resembles a P-Cygni profile, similar to that recently 
found by Nardini et al. (2015) in PDS456. In the latter case the 
emission line peaks however at ~ 7 keV rest frame, suggest¬ 
ing highly ionized Fe emission, unlike the case of Mrk 231. The 
deficit of counts around 7 keV is also very well visible in Figures 
4 and 5 of T14. We replaced the narrow emission line with a rel¬ 
ativistic disk-line with inner and outer radii fixed at 10 and 1000 
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Fig. 13. Left panel: integral radial profiles (histogram) and profiles in annuli (points) of the outflow rate M 0 f (left axis), and v max 
(right axis). Middle panel: radial profile of Ekin,OF- Right panel: radial profile of the outflow/molecular disk mass, MoF/Mdisk- 


gravitational radii, and spectral index of the power emissivity 
law fixed at -2. The^ 2 improves to 470.0 for 467 dof. The equiv¬ 
alent width of the disk-line is 270 eV for the Chandra spectrum 
and around 130-140 eV for the four NuSTAR spectra. Residuals 
still show a deficit of counts around 7 keV. 


We then added to the model a Gaussian absorption line, to 
account for this deficit of counts. The x 2 further improves to 
452.45 for 464 dof when the normalization of the gaussian ab¬ 
sorption line is kept linked in the five spectra. The A^ 2 between 
the model without and with the gaussian absorption line is 17.55 
for three additional dof and no systematic residuals are seen in 
the iron line complex. Figure 15 shows the^ 2 contours for the 
absorption line normalization and energy. The detection of the 
absorption line is significant at better than the 99.9% confidence 
level (or 3.5 <t). To confirm the statistical significance of the ab¬ 
sorption feature we run a Markov Chain Monte Carlo simulation 
by using the XSPEC chain command. This generates a chain 
of set of parameter values describing the parameter probability 
distribution. The chain is converted into probability distribution 
using the XSPEC margin command. 


The Ax 2 between the models with and without the absorp¬ 
tion line for the single Chandra spectrum is 11.6, while it is 7.3 
between the two models for the four NuStar spectra. The en¬ 
ergy of the line is also similar in the two instruments, 7.1 ± 0.5 
keV for the Chandra spectrum and 7.1^4 keV for the NuSTAR 
spectra. We conclude that the line is present with similar statis¬ 
tical significance in both instruments. The equivalent width of 
the absorption line is about twice in the NuSTAR spectra than 
in the Chandra spectrum. Fitting simultaneously the Chandra 
and NuSTAR spectra, but allowing the normalization of the ab¬ 
sorption line to be different in the Chandra and NuStar spectra, 
reduces the^ 2 to 448.24, i.e. the A^ 2 is only 4.2. For the sake of 
simplicity, in the following analysis we keep the normalization 
of the absorption line linked between the Chandra and NuSTAR 
spectra. 


To further test the robustness of the detection of the absorp¬ 
tion line at about 7.1 keV, we performed several additional spec¬ 
tral fits. First, we limited the band used to fit the NuSTAR data 
to 30-40 keV, up to which the source is detected with a signal to 
noise better than 3cr in bins smaller than 10 keV. The results were 
perfectly consistent with those obtained fitting the NuSTAR data 
in the full band, up to 79 keV. Second, we considered each single 
Chandra spectra (eight different observations) and fitted them 
simultaneously. The results were consistent with those obtained 



Energy (keV) 


Fig. 14. The ratio between data the best model including two 
thermal components, a power law component and a narrow gaus¬ 
sian emission line at 6.18 keV (6.44 keV rest frame, iron K a). 
Black triangles = Chandra data; red dots = NuSTAR data. A 
strong excess around 6 keV and a deficit of counts around 7 keV 
are seen. Note that the NuSTAR data, even taking into account 
the lower spectral resolution, follow the Chandra data. 


fitting the single Chandra spectrum obtained by adding together 
the eight observations. Third, we fit the Chandra and NuSTAR 
data with a model including an absorption edge at energy >7.1 
keV, rest frame. The x 2 obtained linking the optical depth r of 
the edge between the Chandra and NuSTAR data is 470.5 for 
465 dof, and the resulting r is consistent with zero. Allowing 
the optical depth of the edge to be different in the Chandra and 
NuSTAR spectra produces a x 2 =45 8.3 for 464 dof, again much 
higher than the^ 2 = 448.24 obtained using an absorption line. 

We fitted the Chandra and NuSTAR spectra with a model 
including two thermal components (T1 and T2), a power law 
nuclear component (PL) and a disk line component. The latter 
two components are modified by two ionized absorbers (a low 
ionization absorber to account for the cold-warm gas along the 
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line of sight, and a highly ionized absorber to account for the 
fast wind). The x 2 is again good, 452.4 for 463 dof. The best 
fit parameters for all components are given in Table [5] Figure 


15 shows the^ 2 contours for the highly ionized absorber N h. 


ionization parameter and redshift. From the absorber redshift we 
derive a velocity of the absorbing gas of Vufo = -20000^^ 


km s -1 . The covering factor of the nuclear wind is consistent 
with being 1. The statistics only allows to obtain a lower limit of 
0.7 at the 67% confidence level, as derived from the partial cov¬ 
ering model (e.g. from the ratio of absorbed/unabsorbed power 
law, Table [5]). In principle, the wind covering factor may be de¬ 
termined by a detailed analysis of the P-Cygni profile of the iron 
line feature, if the emission and the absorption come from the 
same expanding envelope (as in Nardini et al. 2015). Here we do 
not attempt such a detailed analysis, because the statistics does 
not allow us to distinguish between line emission from a wind 
and line emission from a relativistic accretion disk. We limit 
ourselves to note that a prominent P-Cygni profile (e.g. strong 
emission in addition to strong absorption, Figure [14]) points to¬ 
ward a large wind covering factor, provided that both features 
are produced by the same gas. 


The comparison between our best fit model and the preferred 
model of T14 is not straightforward because we used a simpler 
ionized/partial covering absorber model for the nuclear obscu¬ 
ration, while they adopted the MYTORUS model. Furthermore, 
we used a two temperature MEKAL model for the starburst com¬ 
ponent, while T14 use two temperature MEKAL model plus a 
power law component and a 6.7 keV Fe line to account for con¬ 
tribution of a population of nuclear High Mass X-ray Binaries 
(HMXB). Finally, we used a disk Fe K a line, while this com¬ 
ponent is automatically included in the MYTORUS model in 
T14. The best fit spectral index and luminosity of the nuclear 
X-ray component are similar in the two fits, T = 1.47 and 
L(2-10keV) = 3xl0 42 ergs/s versus T = 1.4andL(2-10keV) ~ 
4 x 10 42 ergs/s. Our parametrization of the starburst component 
provides a warm component (T=0.89 keV, similar to T14), plus a 
hot component (5.5 keV, which naturally accounts for the Fe 6.7 
keV line). We tried to add a flat (r = 1.1) power law component 
to account for a HMXB population, but this does not improve 
significantly the^ 2 , and its normalization is < 30 times that of 
the AGN component (similarly to T14). We fit a Compton thin 
absorber with column density ~ 4 x 10 22 cm -2 , very low ioniza¬ 
tion parameter, uniformly covering the X-ray source, while T14 
fit the data with two neutral absorber of similar total column den¬ 
sities. In conclusion, despite significant differences in the details 
of the adopted models, the broad band X-ray spectral reconstruc¬ 
tion is quite similar, and characterized by a) a main AGN power 
law component with a rather flat spectral index; b) a Compton 
thin, moderately ionized absorber of ~ 4 x 10 22 cm -2 , uniformly 
covering the source; c) a warm thermal component and a hot 
component (or a flat power law plus a 6.7 Fe line) to account 
for the nuclear starburst. The detection of the absorption feature 
at about 7 keV is clearly not an artifact of the different spectral 
modeling, see discussion above, and it is well visible also in the 
residuals of Fig. 3 of T14. 


The other outflow parameters can be computed from the best 
fit parameters as follows. The minimum distance of the outflow¬ 
ing gas can be estimated from the radius at which the observed 
velocity equals the escape velocity, r m i n = 2 GM B h/Vu F 0 = 
5 x 10 15 cm for a black hole mass of 8.7 x 10 7 M© (Kawakatu 
et al. 2007) and the lowest value of the outflow velocity. The 
maximum distance r max of the outflowing gas can be computed 
from the definition of the ionization parameter f , 


N H(/nin)^(min) 


= 3.0 x 10 16 cm 


(4) 


where L ion ~ 10 43 ergs/s is the luminosity integrated in the 13.6 
eV - 13.6 keV range. The mass outflow rate at r min and r max can 
then be estimated using: 


M ufo = 0.%nm p N u H FO v UFO rf g = 0.3 - 2.1 M 0 yr~ l (5) 

where N^ F0 = 2.7 x 10 23 cm -2 , and f g is a geometrical factor 
introduced by Krongold et al. (2007) and statistically estimated 
~ 2 by Tombesi et al. (2010). Should the mass outflow rate be 
significantly larger than the value given above, a much deeper 
absorption line should be visible in the X-ray spectra (for reason¬ 
able values of f and v out ). From the mass outflow rate and veloc¬ 
ity we can estimate the momentum flux, Pufo = (0.4-2.7) x 10 35 
gr cm/s 2 , corresponding to a Ekin,UFO = 3.8 x 10 43 - 2.7 x 10 44 
erg s -1 . This can be compared with the radiation momentum flux 
P ra d = Lboi/c. The bolometric luminosity estimated by Lonsdale 
et al. (2003) and Farrah et al. (2003) is 5 x 10 45 erg s -1 . The mo¬ 
mentum load is therefore Pufo /(Prad) = 0.2-1.6. Recently, Teng 
et al. (2014) measured an intrinsic X-ray luminosity in the range 
0.5-30 keV of ~ 10 43 erg s -1 , about 10 times smaller than previ¬ 
ous estimates. This would either imply an unusually large X-ray 
bolometric correction (« 1000) and hence a very large momen¬ 
tum load, or an overestimated AGN bolometric luminosity based 
on mid-infrared data. 

The detection of an UFO offers the unprecedented opportu¬ 
nity to directly compute the momentum load of the large scale, 
spatially resolved outflow with respect to that of the nuclear 
wind. The Pof of the CO(2-l) outflow depends on velocity and 
is in the range 1.6-4x 10 36 erg cm -2 for v=500 and 800 km s -1 , 
respectively. Figure [16] left panel, shows the momentum load, 
Pof/Pufo, versus the velocity of the molecular and atomic large 
scale outflows detected in Mrk 231 (references for the data are 
in the caption of the figure). The values of Pof/Pufo range from 
« 20 to « 100, i.e. strongly suggesting that the outflow does not 
conserve the momentum, if the UFO is identified with the inner 
semi-relativistic wind pushing the outer molecular outflow (Lapi 
et al. 2005, Zubovas & King 2012, Faucher-Giguere & Quataert 
2012). Figure [l6] right panel, shows the P outflow /Prad versus ve¬ 
locity of the oWnows of Mrk 231 (red symbols). Black symbols 
are the results of Tombesi et al. (2015) for IRAS Fill 19+3257, 
a galaxy exhibiting both an X-ray UFO and an OH outflow, as re¬ 
vealed by Herschel spectroscopy (Veilleux et al. 2013). The solid 
lines represent the expectations for energy-conserving outflows 
with nuclear wind velocity vufo ± lo - and with covering factor 
/ = 1 for both the nuclear winds and the molecular outflows. 

5. Discussion and conclusions 

(i) Molecular disk. We present the best spatial resolution and 
sensitivity CO map of the molecular disk of Mrk 231 obtained 
so far. In addition to the previously known main rotation of the 
disk, which occurs along nearly an east-west direction, we find 
evidence for a tilted inner nuclear disk on scales of 100 pc, more 
inclined on the line of sight. The CO gas distribution can be mod¬ 
eled by an outer exponential ring with radius 400 pc, PA = -12 
deg, and inclination i = 36+10 deg, plus an inner exponen¬ 
tial disk with outer radius 200 pc, PA = 84 deg, and i = 58 
deg on the line of sight. The inner, warped component matches 
the size, PA and inclination of the vibrationally excited HCN(3- 
2) v=l,f (Aalto et al. 2014). The parameters of the inner tilted 
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Fig. 15. Left panel: 67%, 90% and 99% iso-^ 2 confidence levels derived from the fitting of the merged X-ray spectra for the 
absorption line normalization and energy (left panel), for the column density, Nh, and redshift of the highly ionized absorber 
(middle panel), and for N H and ionization parameter, of the highly ionized absorber (right panel). 
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Fig. 16. Left panel: the momentum boost Pof/Pufo for the known molecular and atomic outflows of Mrk 231. References are Aalto 
et al. (2014) for HCN, Gonzalez-Alfonso et al. (2014) for OH, Rupke & Veilleux (2011) for NalD, this work for CO. Right panel: 
the ratio of Pof/P rad versus velocity of the outflows of Mrk 231 (diamond: this work, stars: other works, see references above). 
Filled squares show the measurements for IRAS F11119+3257 (Tombesi et al. 2015). Red and black solid lines are the expectations 

for energy conserving outflows, = y 24 , with v w i n d = vufo ± 1 cr for Mrk 231 and IRAS FI 1119+3257, respectively. 


disk are also consistent with the OH mega-maser emission seen 
with the VLBI, which traces a disk with radius 30 to 200 pc, 
PA 56 deg and inclination 56 deg (Klockner et al. 2003). We 
can hardly estimate the mass of the inner tilted CO disk, but 
this detection shows that in addition to very dense gas traced by 
HCN(3-2) V 2 = 1/, gas with density of the order 10 5 cm -3 is 
also present in the very inner regions. The existence of a warped 
inner disk with radius 170-250 pc was suggested also by Davies 
et al. (2004) based on stellar absorption studies. They pointed 
out that the warp should first occur in the gaseous phase and so 
be transferred to the stars formed in situ. 

(ii) Molecular outflow. The CO(2-l) high angular resolution 
data indicate that the molecular outflow has a size of at least 
1 kpc. The bulk of both receding and approaching outflowing 
gas are located within -400 pc from the nucleus, and peak ~0.2 
arcsec south-west off the nucleus. Extended, redshifted emission 
with lower surface brightness is seen north-east off the nucleus 
out to - 1 kpc (see Fig. 3). 

The outflow is seen in all examined directions around the 
AGN (see Fig. 6). It is more prominent, however, along the 


south-west to north-east direction, suggesting a wide-angle 
likely biconical geometry. The approaching gas shows a gradient 
in projected velocity, with larger negative speed at the position 
of the AGN and decreasing velocity south-west off the nucleus. 
This allows to exclude a geometry where the approaching out¬ 
flow expands with constant speed within an uniformly filled cone 
with axis located in the plane of the sky (as in this case no speed 
gradient would be observed) nor with axis parallel to our line of 
sight (in this case the maximum projected speed would be cen¬ 
tered on the nucleus and decrease outwards symmetrically). An 
outflow expanding with uniform velocity between the observer 
and the molecular disc, in a cone with axis inclined with respect 
to our line of sight, would instead produce a velocity gradient 
as the observed one for the approaching gas, but it would also 
imply a similar gradient for the receding part, which instead is 
not observed. The projected velocities can be probably qualita¬ 
tively be explained by an outflow with conical or lobe-like ge¬ 
ometry with a slightly larger receding component (based on Fig. 
2 and 3). However a firm and consistent picture explaining both 
velocity fields seen in Fig. 9 cannot be simply ascribed to orien- 
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Fig. 17. 3D cartoon of the molecular disk and outflow geometry in Mrk 231. The line of sight to the observer is perpendicular to 
the page. Red and blue colours represent the receding and approaching material for each component, respectively. The molecular 
disk consists of an inner disk, tilted by 58 deg on the line of sight, with a scale radius of 200 pc, plus and outer component with 
inclination 36 deg and scale radius 400 pc, with the line of nodes oriented in the east-west direction approximately. The molecular 
outflow is rendered by a wide-angle biconical-like component plus an isotropic one, represented here as a sphere for displaying 
purposes. For a detailed description of the cartoon see Section 5. 


Table 5. Best fit model to the Chandra and NuSTAR merged X 
ray spectra 

Parameter 

Best fit value 


Units 

KT(1) 

0.89±0.05 


keV 

Normalization 1) 

(4.9+0.5) x 10" 6 

P h 

cm -2 s -1 keV -1 

KT(2) 

5.5±0.7 


keV 

Normalization^) 

(5.7 ± 2.1) x 10 -5 

P h 

cm -2 s -1 keV -1 

Nff(l) 

(4.2+0.3) x 10 22 


_2 

cm 

Covering Factor(l) 

0 Q7+0.03 
^’ y ' -0.04 



log?(l) 

-0 54+0- 16 

-0.04 


erg cm s -1 

N*(2) 

(2.7+fJ) x 10 23 


cm" 2 

Covering Factor(2) 

PO-0.3 



log?(2) 

3 55 +0 - 2 


erg cm s _1 

F(PL) 

1.47+0.1 



Normalization(PL) 

(1.9+0.2) x 10 -4 

P h 

cm -2 s _1 keV -1 

Disk Line Energy 

6.05+0.04 


keV 

Disk Inclination 

35+5 


deg 

Normalization(Line) 

(3.0+0.5) x 10" 6 

P h 

cm -2 s 1 keV -1 

AT 2 (dof) 

452.4 (463) 




tation effects and probably inhomogeneities and asymmetries in 
the outflow can account for them. 

We also suggest that the redshifted emission seen north-east 
off the nucleus is tracing the receding part of the cone, which is 
located behind the molecular disk with respect to the observer. 
Interestingly, Krabbe et al. (1997) detected at approximately the 
same location a narrow Pa a emission approaching with a speed 
of 1400 km s -1 with respect to the systemic one. The molecular 


outflow centroids, traced by the red and blueshifted wings of 
CO(2-l), are consistent with those seen in HCN (Aalto et al. 
2014). 

We find that Mof decreases with increasing distance from 
the nucleus in all directions (see Fig. 11), while the radially in¬ 
tegrated v max stays roughly constant out to > 1 kpc. This im¬ 
plies, based on equation [3] that either the outflow average den¬ 
sity, pof> or the outflow filling factor decrease from the nucleus 
outwards approximately as r~ 2 (or slightly steeper). This sug¬ 
gests different possible scenarios: (i) a large part of the gas with 
velocity of 300-600 km s -1 leaves the flow during its expan¬ 
sion. If the molecular clouds are pressure-confined, they will dis¬ 
solve out in the wind, and CO may be efficiently photo dissoci¬ 
ated by the UV radiation, since self shielding will be strongly 
reduced at low densities. Atomic carbon, either in neutral or 
ionized phase, can then, in principle, continue expanding out 
to larger distances at a constant velocity, as it is observed in 
SDSSJ1148+5251, where most of the fast [CII] is located far 
from the AGN, and extended on scales of a few tens kpc (Cicone 
et al. 2015). Intriguingly, broad wings extending up to ±1000 
km/s can be seen in both [CII] 158 pm and [Nil] 205 pm emis¬ 
sion lines observed by Herschel in Mrk 231 (Fischer et al. 2010). 
Unfortunately, Herschel PACS does not provide sufficient angu¬ 
lar resolution to assess whether the fast [CII] and [Nil] gas are 
more extended than CO. A promising alternative would be to 
use the forbidden 3P fine structure line of atomic carbon ([Cl] 
line at rest frame frequency 492.161 GHz), which is accessible 
by ALMA, on southern analogs of Mrk 231. (ii) The opening 
angle subtended by the paths of least resistance is relatively nar¬ 
row, or (iii) only the highest velocity gas has been able to reach 
distances > 1 kpc. The latter would imply an age of the outflow 
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of - 1 Myr (assuming a constant expansion), which is 1-10% of 
a typical AGN duty cycle. The rate of kinetic energy carried by 
the outflow is — 1-2% of the AGN bolometric luminosity, and 
stays approximately constant out to 1 kpc. The decrease of the 
local Ekin,oF at larger radii can be explained as a consequence of 
decreasing M 0 f and constant v max . 

A 3D rendering of the geometry of the molecular disk and 
outflow is shown in Figure [IT] The molecular disk consists of an 
inner tilted disk plus and outer component with approximately 
an east-west orientation. The observed direction of the rotation 
for these two rotating components is shown by the colour scale 
(red for the receding, blue for the approaching material). The 
molecular outflow is sketched out by a wide-angle biconical-like 
component plus a spherical one. The two cone-like components 
have been placed in order to match the south-west to north-east 
elongation seen in the CO(2-l) maps. The southern cone is lo¬ 
cated between the observer and the disk, inclined with respect to 
the line of sight, and shows an approaching and a brighter reced¬ 
ing gas emission. The northern cone is located behind the molec¬ 
ular disc, and it is mainly receding, as testified by the redshifted 
velocities observed at this position. Since we cannot draw firm 
conclusions on the expansion velocity along the plane of the sky, 
we mention the possibility that the isotropic outflowing compo¬ 
nent can also be explained by receding and approaching outflows 
along the line of sight. It is worth noting that the neutral gas out¬ 
flow reported by Rupke & Veilleux (2011) extends over a much 
larger region (-3 kpc) than that of the molecular gas traced by 
CO. 

(iii) Molecular outflow and disk global energetics. A com¬ 
plementary way to assess whether the outflow has a signifi¬ 
cant impact on the molecular disk is comparing their energet¬ 
ics, keeping in mind, however, that these energy estimates are 
affected by large uncertainties. The total energy of the molecular 
disk (quiescent component) is given by the sum of the gravita¬ 
tional potential energy, the rotational energy, the turbulent en¬ 
ergy and the thermal one. The gravitational potential energy is 
Egrav = GM dyn M gas /R, where M dyn and M gas are the dynamical 
and gas masses, G is the gravitational constant and R is the radius 
of the molecular disk (400 pc from our fit). The gas mass within 
this radius is M gas = 1.8 x 10 9 M©. The dynamical mass, mod¬ 
ulus the inclination of the disk, is M^ yn sin 2 (i) = 5.5 x 10 8 M 0 , 
which gives E grav sin 2 (i) = 1.85 x 10 56 erg. 

The rotational energy of the molecular disk is E rot = 
\M gas v 2 = 2.1 x 10 57 erg, where v = 345 km s -1 is the rotational 
velocity of the disk. The contribution of the turbulent energy is 
small and can be expressed as: E turb = \M gas dv 2 = 2.9 x 10 55 
erg, where dv is the turbulent velocity dispersion in the disk 
(which is likely of the order of ~ 40 km s -1 ). The thermal energy 
can be expressed as E t h e rm = n k T V = 1.2 x 10 52 erg, assum¬ 
ing a temperature of the gas of T = 60 K based on CO(l-O) 
observations, an average gas density n = 3600 cm -3 , and a disk 
with radius 400 pc and thickness 23 pc (Downes & Solomon 
1998). The total energy budget of the molecular disk is therefore 
E d isk = 2.6-8.3x10 57 erg, adopting a disk inclination of 36 and 
10 deg, respectively. 

The kinetic energy of the outflow is given by its bulk motion. 
Based on the mass of the CO(2-l) outflow, we find Ekm,oF = 
2.35 x 10 57 erg, adopting an average velocity of 800 km s -1 and 
an outflow mass of 3.7 x 10 8 M 0 The turbulent contribution is 
negligible if we assume dv turb , 0 F = dv turb 4 isk ~ 40 km/s. On 
the other hand if in the outflowing molecular clouds dv turb ~ 
250 km/s (e.g. Williamson et al. 2014), the turbulent energy can 
contribute an additional ~ 10% to the outflow energy budget. 


Given all the uncertainties in these estimates, we conclude that 
the total energy of the outflow Eof is of the same order as E d i S k. 

(iv) UFO and molecular outflow energy and momentum. We 
find a ratio = 0.4 - 3.8, suggesting that most 

of the nuclear wind kinetic energy is transferred to mechanical 
energy of the kpc scale outflow, which is thus undergoing an 
energy conserving expansion, as predicted by the most popular 
theoretical models (Faucher-Giguere & Quataert 2012, Zubovas 
& King 2012). In particular, the prediction of models for energy- 
conserving outflows is Ekin,uFO =/xEkin,OF, where / is the frac¬ 
tion the nuclear wind energy that is deposited to the large scale 
outflow, and ranges from 0.5 to 1 in the models of Faucher- 
Giguere & Quataert (2012) and Zubovas & King (2012), respec¬ 
tively. We stress that our work does not make any assumption on 
/. In fact, we measure 

Pqf _ vufo 

Pufo v of 

from which we derive / - 1 (see also Figure [T6|). The mass 
loading factor, defined as // = sJ{Mof I Mjjfo ) in Zubovas & 
King (2012), is « 20 at the outer boundary of the molecu¬ 
lar outflow (1 kpc), which also supports the two stage accel¬ 
eration scenario. We find that Ekm,UFo/Lb 0 i,AGN = 1-5% and 
Ekin,OF/Lboi,AGN = 1 - 3%, in agreement with the requirements 
typically assumed by the most popular feedback models (Lapi et 
al. 2005, Menci et al. 2008, Di Matteo 2005, Gaspari et al. 2011, 
2012 ). 

The momentum load of the nuclear wind, 
Pufo/(L boi ,agn /c) = 0.2 - 1.6, agrees with the predictions 
for radiatively accelerated winds with scattering optical depth 
~ 1 (King & Pounds 2003). These observations offered the first 
opportunity to compare the momentum boost of a spatially re¬ 
solved outflow with respect to the nuclear wind. Specifically, we 
find Pof/Pufo = [16 - 113] at 1 kpc based on CO(2-l), and in 
the range 8 -100 for the outflows traced by HCN, OH and NalD, 
although with large uncertainties (Figure [16]). Accordingly, this 
result supports the two phase outflow acceleration mechanism 
scenario, whereby the momentum of the large scale outflow is 
boosted compared to that of the nuclear semi-relativistic wind. 
We remark that the nucleus of Mrk 231 is radiating close to the 
Eddington limit (L bo i,AGN/LEdd = 0-46 for a black hole mass of 
8.7 x 10 7 M 0 ), matching again the requirements of models for 
driving massive outflows. 

Tombesi et al. (2015) recently found a similar result for the 
galaxy IRAS FI 1119+3257, which hosts both an UFO and a 
massive molecular outflow detected through the OH absorption 
line by Herschel (Veilleux et al. 2013). They derived PoF/Prad ~ 
12, estimating / = 0.2. It is worth noting that both Tombesi et 
al. (2015) and this work strengthen each other, supporting the 
AGN-driven energy-conserving outflows scenario and providing 
constraints to the models of AGN feedback. 

Sensitive and high angular resolution observations with 
ALMA and NOEMA are needed to further constrain the physics 
of quasar-driven outflows, and their impact in galaxy trans¬ 
formation. In particular, constraining any large-scale molecu¬ 
lar outflow in quasars with well-studied UFO (e.g. PDS 456, 
PG 1211 + 143, Nardini et al. 2015, Pounds 2014) is a compelling 
experiment to be pursued with ALMA. On the other hand, under¬ 
standing the impact of outflows on the cosmological evolution of 
galaxies requires a systematic approach, i.e. a search for nuclear 
molecular outflows out to z ~ 2 in unbiased, mass-selected sam¬ 
ples. 
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